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INTRODUCTION 


This  report  documents  modifications  which  were  incorporated  in  the  EPIC- 3 
[5,6]  computer  program  (Elastic- Plastic  Impact  Computations  in  Three 
Dimensions)  in  order  to  treat  impact  and  penetration  with  erosion.  EPIC- 3  is 
a  finite  element  program  with  explicit  time  integration  which  is  primarily 
intended  for  simulation  of  solids  subjected  to  short  intense  loads,  such  as  in 
impact  or  explosive  detonations.  A  variety  of  material  laws  including 
elastic-plastic  solids,  concrete/geological  materials  and  explosives  are 
included. 

As  part  of  t*:is  modification,  two  major  features  were  added: 

1 .  a  hexahedral  element  with  one-point  quadrature  and  hourglass  control; 

2.  an  algorithm  for  treating  projectile- target  interaction  in  situations 
where  material  erosion  can  occur  arbitrarily  in  the  target  or  projectile. 

The  hexahedral  finite  element  is  an  isoparametric  element  with  8  nodes  and 
24  degrees  of  freedom.  The  stresses  and  strains  are  evaluated  only  at  a 
single  point  in  the  element,  which  makes  certain  physically  unrealistic  modes 
of  deformation,  known  as  hourglass  modes,  possible.  These  are  avoided  here  by 
the  use  of  hourglass  control. 

A  key  attribute  of  the  interaction  algorithm  is  that  it  requires  no 
definition  or  tracking  of  sliding  interfaces.  Instead,  the  interaction  is 
handled  by  operations  on  slave  nodes  and  master  elements.  Because  of  this 
feature  of  the  algorithm,  the  erosion  of  an  element  requires  no  redefinition 
of  the  interface  and  thus  avoids  the  complexity  associated  with  sliding 
interfaces  in  these  situations.  The  hexahedral  element  was  incorporated 
primarily  because  it  simplifies  the  new  interaction  algorithm.  However,  it 
also  increases  the  speed  of  the  computer  program  and  avoids  the  excessive 
stiffness  of  the  tetrahedra. 

The  element  is  described  in  Section  2,  the  interaction  algorithm  in 
Section  3.  Section  3  describes  some  sample  problems  which  were  solved  by  this 
algorithm.  The  modified  input  manual  is  given  in  Section  4.  Equations  in 
this  report  are  numbered  by  section  and  sub- section,  so  that  each  equation 
number  consists  of  3  numbers;  the  section  number,  the  sub-section  number,  and 
the  number  of  the  equation  with  the  sub- section.  However,  when  referring  to 
an  equation,  the  first  number,  the  section  number,  is  included  only  when  the 
equation  is  in  a  different  section. 


Section  1 


HEXAHEDRAL  ELEMENT 


1  .1  Overview  and  Notation 


The  hexahedral  element  consists  of  8  nodes  and  6  sides  as  shown  in 
Fig.  1  .  The  positions  of  the  nodes  are  completely  arbitrary  although  the 
local  node  numbering  must  be  in  accord  with  the  following  convention: 

i)  nodes  1  to  4  must  be  placed  sequentially  on  a  single  surface  so  that  the 
thumb  of  the  right  hand  points  to  the  interior  of  the  element  when  the 
fingers  follow  the  nodes  in  sequence; 

ii)  nodes  5  to  8  must  each  be  connected  by  edges  to  nodes  1  to  4,  respec¬ 
tively;  i.e.,  nodes  1  and  5,  2  and  6,  3  and  7  and  4  and  8  must  each  lie  on 
the  same  edge . 

The  element  is  an  isoparametric  element  with  the  displacement  field 
defined  in  terms  of  the  coordinates  of  reference  cube  in  £,  n,  C  space,  see 
Fig.  1 .  In  terms  of  the  reference  coordinate  system,  the  displacement  and 
velocity  field  are  trilinear,  with  none  of  the  coordinates  appearing  in  the 
polynomial  in  a  power  higher  than  linear.  The  resulting  strain  fields  are 
then  bilinear. 

However,  only  one  quadrature  point  is  used  in  evaluating  the  strain  and 
stress  fields  in  the  element.  This  implies  that  for  purposes  of  evaluating 
the  nodal  forces,  the  strain- rates  and  stresses  are  considered  constant. 

The  assumption  of  a  constant  stress  field  implies  that  certain  deforma¬ 
tion  modes  of  the  element  will  not  be  resisted  by  nodal  forces;  this  pheno¬ 
menon  is  known  as  hourglassing;  see  Ref.  (1].  To  avoid  the  severe  mesh 
distortions  brought  about  by  hourglassing,  an  hourglass  procedure  developed  by 
Flanagan  and  Belytschko  [2]  will  be  used.  This  hourglassing  procedure  is 
advantageous  in  that  it  does  not  compromise  the  formal  consistency  of  the 
resulting  difference  equations,  so  convergence  is  not  impaired;  see  Ref.  [1]. 

In  the  following,  we  first  give  the  fundamental  equations  for  the  hexa¬ 
hedral  element.  Two  types  of  £  matrices  for  one- point  quadrature  have  been 
incorporated;  one  is  based  on  centroidal  evaluations  of  the  B  matrix  as  in 
Ref.  [3]  the  second  is  based  on  the  uniform  strain  procedure  given  in 
Ref.  [21.  The  centroidal  method  is  cheaper  but  loses  accuracy  as  the  element 
distorts,  in  particular,  if  the  element  is  degenerated  to  a  pentahedron  or 
tetrahedron,  it  becomes  very  inaccurate,  whereas  the  Flanagan- Belytschko 
formulas  [2]  remain  exact. 

After  the  development  of  the  ^-matrices,  the  hourglass  control  procedure 
is  given.  Then  the  formulas  for  estimating  the  stable  time  steps  based,  on 
Ref.  [4],  which  are  used  for  the  hexahedron,  are  presented. 


Standard  indicial  notation  will  be  used  in  this  report.  Lower  case  Latin 
indices  pertain  to  components  and  have  a  range  of  three;  when  they  are 
repeated  a  summation  over  their  range  is  implied.  Upper  case  indices  pertain 
to  nodes  of  an  element  and  have  a  range  of  8.  The  range  of  Greek  letters  is 
specified  whenever  they  are  used.  Commas  denote  derivatives  with  respect  to 
the  spatial  variables,  i.e.  f,^  *  Sf/Sx^. 

Where  convenient,  matrix  notation  is  used.  Matrices  are  denoted  by  tilde 
subscripts,  such  as  p.  Lower  case  letters  designate  column  matrices  or 
vectors,  whereas  upper  case  letters  designate  rectangular  matrices.  The 
superscript  "T"  denotes  a  transpose,  as  in  n  . 


1.2  Mapping,  Displacement  Field  and  Fundamental  Equations 


The  mapping  between  the  physical  and  reference  domains  of  the  element, 
which  are  shown  in  Fig.  1 ,  are  given  by 

xi  a  N  ^(  £ ,  n,  £)  (1.2.1) 

where  £,  n,  £  are  the  coordinates  in  the  reference  cube  of  dimension  2x2x2 
with  the  origin  at  the  centroid  of  the  element. 

The  shape  functions  Nj  are  given  by 

Njtc,  n,  c)  -  J  (1  +  £j£)  (1  +  nIn)  (1  +  CICI>  (1.2.2) 

where  £  ,  n,  £  are  the  coordinates  of  node  I  in  the  reference  domain,  which 
are  all±  1 . 


The  displacement  and  velocity  fields  in  the  element  are  given  by  the  same 
shape  functions  since  the  element  is  isoparametric,  hence 

u^  =  u^  Nj  (1  .2.3a) 

v^  ■  v^j  Nj  (1 .2.3b) 

where  u^  and  v^  are  the  displacement  and  velocity  fields  and  u^  and  v 

the  components  of  the  nodal  displacements  and  velocities  of  node  I. 


The  velocity  strains  are  given  by 


e 


ij 


aoi 

3x.  Vil 
3 


(1.2.4) 


and  the  spin  rates  by 


w .  *  e .  . ,  W ,, 

i  ijk  jk 


(1  .2.5a) 
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Fig.  1  Element  reference  and  physical  domain 


(1  .2.5b) 
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The  element  nodal  forces  are  given  by  £2] 
3h 

fil  "  /  3x7  °ij  dV 

„e  3 


(1 .2.6) 


where  a . .  are  the  stresses 
of  shape3 function  matrix, 
by 


.  If  we  let  B  be  an  average  value  of  the  gradient 
then  in  one- point  quadrature  Eq.  (2.6)  is  replaced 


£il  ■  (Bjl  "ij  *  \l  5Ol)  • 

Only  the  first  term  on  the  right  side  of  (2.7)  originates  from  (2.6); 

Qa^,  a  =  1  to  4,  are  the  generalized  antihourglass  stresses  and  y  are  the 
hourglass  operators  which  are  defined  later  and  must  be  added  to  control 
hourglass  modes. 


1  .3  B  -  Matrix 

Two  forms  of  the  B- matrix  are  incorporated  in  the  program.  In  the  form 
given  in  Ref.  [2]  , 

B11  “  12VB11  3  y2(z63  +  z54}  +  y3Z24  +  y4{*38  +  Z25) 

+  y5(z86  +  z42 }  +  y6Z52  +  y8Z45  (1.3.1a) 

where 

XIJ“XI-XJ'  yU-yi-yj'  ZIJ“ZI-ZJ*  (1.3.1b) 


The  other  terms  of  the  B  matrix  in  row  1  are  obtained  by  simply  permuting  the 
nodal  coordinates  according  to  Table  1 . 
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where 


*  72  [(2ys 

-  V 

Z42 

+ 

y2 

(Z53 

+ 

Z54) 

+  y4 

(-Z53 

-  *52)] 

-72  [<y4- 

2y5> 

Z31 

+ 

y3 

(Z54 

+ 

*51  } 

+  yi 

('Z54 

-  *53>] 

-  72  [(yi  - 

2ys> 

Z42 

+ 

y4 

(Z51 

+ 

Z52> 

+  y2 

(-Z51 

"  Z54>} 

*71  f  (2y5 

-  y2} 

Z31 

+ 

yi 

(Z52 

+ 

Z53) 

+  y3 

(‘Z52 

-  *51)] 

-  71  [  (y54_ 

y52): 

Z3l“ 

(y53 

-  y51 

)z 

:42  ” 

(Z51 

Z53 

)y42  +  (Z52 

(2.3.5) 


In  using  these  formulas,  nodes  1  to  4  must  define  a  side  of  an  element  and 
must  be  numbered  so  that  they  are  counterclockwise  when  viewed  from  a  point 
inside  the  element. 


As  soon  as  any  volume  is  found  to  be  negative,  the  slave  node  can  be 
considered  to  definite)y  lie  outside  the  element.  All  6  pentahedral  volumes 
must  be  positive  if  the  nodes  are  inside  the  element. 


2.4  Normal  Directions 


An  important  ingredient  in  defining  the  interaction  of  the  slave  nodes 
and  master  elements  is  that  any  transfer  of  momentum  which  occurs  between  the 
target  and  penetrator  (other  than  that  due  to  friction,  which  is  not  con¬ 
sidered  here)  should  be  in  directions  normal  to  the  interface.  For  this  pur¬ 
pose,  the  normal  vectors  must  be  available  when  the  positions  of  the  slave 
nodes  are  adjusted.  Since  an  interaction  surface  is  never  defined,  it  is 
necessary  to  construct  normals  in  an  alternative  manner. 

The  assembly  procedure  of  the  finite  element  method  provides  a  very 
natural  and  concise  way  of  computing  these  normal  vectors.  We  will  first 
present  the  procedure  in  a  two  dimensional  setting  so  that  its  ingredients  may 
be  understood  and  visualized  more  readily  and  then  present  the  three 
dimensional  equations. 

The  basic  idea  of  this  method  is  that  the  normals  for  each  side  of  the 
element  are  computed  and  added  component  by  component  into  a  normal  vector 
array  according  to  node  numbers.  The  procedure  is  illustrated  in  Fig.  6. 

Note  that  on  interior  nodes  the  assembled  normal  vectors  essentially  cancel, 
so  their  components  are  very  small  or  zero.  On  exterior  nodes,  the  normal 
vectors  point  out  from  the  domain  with  a  direction  which  reasonably  approxi¬ 
mates  a  normal  to  a  surface  on  the  edge  of  the  domain.  In  particular,  the 
normal  at  a  corner  takes  an  average  direction  of  the  surfaces  that  meet  at  the 
corner . 

In  three  dimensions,  the  procedure  consists  of  the  following:  for  each 
side  of  the  element  with  local  nodes  1,  2,  3,  4,  a  vector  normal  to  the  side 
is  approximately  computed  by 
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ZK  “  2j)‘ 


'  *K  "  ~XJ  '  =  (*K  "  XJ^  +  "  yj3  +  (’ 

and  where  the  correspondence  between  I,  J  and  K  is  given  by 


(2.3.3) 


I  K  J 

1  7  1 

2  8  2 

3  6  4 

4  3  5 


The  square  of  the  radius  is  used  in  all  computations. 

All  slave  nodes  which  are  located  in  cells  which  are  occupied  by  nodes  of 
the  master  element  are  processed  to  check  whether  they  are  in  the  element.  If 
an  element  node  occurs  in  a  cell  which  has  a  leading  zero  in  LOCSLA,  then  the 
cell  contains  no  slave  nodes,  so  no  checks  need  to  be  made.  Otherwise,  the 
process  of  determining  whether  a  slave  node  occurs  in  the  element  begins. 

In  order  to  reduce  computations,  the  slave  nodes  which  are  sufficiently  far 
from  the  master  element  are  first  eliminated  from  consideration  by  the  radius 
check . 


If  a  slave  node  is  within  the  radius  Rg  of  the  master  element,  the  more 
exact  and  time  consuming  checks  are  made  on  the  slave  to  see  whether  it  is 
within  the  element.  If  the  slave  node  is  within  the  element,  the  position  of 
the  slave  node  is  adjusted  and  the  corresponding  momentum  is  transferred  to 
the  appropriate  nodes  of  the  master  element.  In  addition,  the  slave  node 
number  is  deleted  from  the  LOCSLA  list  so  that  this  slave  node  is  not  checked 
against  other  elements.  The  details  of  this  are  given  in  Section  2.5. 

Once  a  slave  node  has  passed  the  simple  checks  for  possible  penetration, 
definitive  check  is  made  as  to  whether  the  node  is  within  the  element.  This 
is  accomplished  by  constructing  six  pentahedra,  each  consisting  of  a  side  of 
the  hexahedron  and  the  slave  node,  as  shown  in  Fig.  5.  If  the  volumes  of  all 
six  pentahedra  are  positive,  the  slave  node  must  be  within  the  element. 

The  volumes  of  the  pentahedra  are  computed  by  using  Eq.  (1.3.4)  with 
nodes  5  to  8  considered  coincident.  This  gives  the  following  formula  for  the 
volume  of  the  pentahedra 


(2.3.4) 


calculations  performed  so  far.  The  number  of  slave  nodes  per  cell  is 
presently  limited  to  40  but  this  can  easily  be  increased  by  changing  the 
dimension  of  the  array  LOCSLA,  as  described  later. 


2.3  Element  Penetration  Detection 


Any  slave  nodes  which  penetrate  an  element  are  treated  during  a  loop 
through  the  elements;  this  loop  is  separate  from  the  element  loop  in  which  the 
nodal  forces  are  updated,  and  considers  only  master  elements.  The  element 
level  steps  consist  of  the  following  for  each  element  e: 

1.  determine  whether  any  slave  node  has  penetrated  element  e; 

2.  if  a  slave  node  has  penetrated  element  e,  place  the  node  outside  the 
element  and  transfer  the  momentum  loss  to  the  element  nodes. 


Note  that  this  procedure  is  carried  out  only  for  master  elements  so  the 
subdomain  of  the  target  that  is  designated  to  consist  of  master  elements  must 
be  sufficiently  large  so  that  no  slave  node  can  penetrate  the  target  without 
penetrating  a  master  element. 

The  algorithm  for  step  1  must  be  streamlined  as  much  as  possible  in  order 
to  insure  that  computer  time  is  not  wasted.  The  cell  scheme  is  used  to  mini¬ 
mize  as  much  as  possible  the  number  of  computations  required.  The  details  of 
the  procedure  are  as  follows. 


The  cell  number  of  each  node  of  the  element  is  determined.  As  shown  in 
Fig.  4,  it  is  possible  for  the  element  to  be  in  several  cells.  (Note  for 
future  work,  it  may  be  possible  to  avoid  this  by  using  an  overlapping  cell 
lattice  so  that  for  any  element  in  the  cell,  all  slave  nodes  contained  in  the 
element  must  be  in  the  same  cell  as  the  centroid  of  the  element.)  The  cell 
numbers  of  the  nodes  are  determined  in  subroutine  LOCNOD,  which  is  the  same 
subroutine  that  is  used  for  determining  the  slave  node  cell  location.  If 
ICELL  *  0,  the  node  is  not  in  the  interaction  zone,  so  it  need  not  be 
considered.  The  cell  numbers  of  all  nodes  of  the  element  are  checked  against 
those  of  the  other  nodes  so  that  no  duplicate  cell  numbers  are  considered  in 
checking  for  slave  node  penetration  of  the  element. 


To  determine  which  slave  nodes  are  in  an  element,  all  slaves  in  the  same 
cell  as  the  element  are  checked.  First  a  rough  check  is  made.  For  this 
purpose,  the  centroid  of  the  element  is  defined  by 
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(2.3.1  ) 


where  x ^  are  the  coordinates  of  node  I.  The  radius  of  the  element  is  defined 
by 


2  i 

R  =  —  max  I  x  —  x,  1  for  1=1  to  4 

e  4  ~K 


(2.3.2) 


Here  and  x  are  the  x-coordinates  of  the  extreme  points  of  the  cell 

domain  and  NCX  is  the  number  of  cells  in  the  x-direction;  the  definitions  of 
the  terms  in  (2.2)  and  (2.3)  are  analogous. 


The  cell  location  of  a  node  with  coordinates  (x,  y,  2)  is  computed  in  two 
steps.  First  the  integer  cell  numbers  in  the  x,  y  and  z  directions  are 
computed  by 


IX  = 

NCX* 

(X  - 

Xmin^ 

/ 

(x 

max 

X  .  ) 
min 

+  1 

(2.2.4a) 

IY  - 

NCY* 

(y  - 

w 

/ 
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max 

y  .  ) 
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(2.2.4b) 

IZ  = 

NCZ* 

(z  - 

z  .  ) 
min 

/ 

(z 

max 

z  .  ) 
min 

+  1 

(2.2.4c) 

These  numbers  are  then  used  to  compute  the  cell  number  of  the  node,  ICELL,  as 
follows.  If  the  node  is  outside  the  cell  domain,  the  cell  number  ICELL  is  set 
to  zero,  i.e. 


ICELL  -  0  if 

IX 

<  1 

or 

IX 

> 

NCX 

or 

IY 

<  1 

or 

IY 

> 

NCY 

(2.2.5) 

or 

IZ 

<  1 

or 

IZ 

> 

NCZ 

Otherwise,  ICELL  is  computed  by 

ICELL-  (IZ-1 )*NCX*NCY  +  (IY-1)*NCX+  IX  (2.2.6) 

The  locations  of  all  slave  nodes  are  stored  in  an  array  K=LOCSLA( IC, J) 
where  IC  is  the  cell  number  and  K  is  the  node  number  of  the  Jth  slave  node  in 
cell  IC.  IF  LOCSLA( IC, 1 )-0,  then  no  slave  nodes  are  located  within  cell  IC. 

The  identification  of  the  cell  locations  of  slave  nodes  is  made  during 
the  time  integration  loop.  As  the  new  positions  of  slave  nodes  are  calculated 
(prior  to  any  adjustment  for  interaction  with  master  elements),  the  cell 
number  is  obtained  by  the  subroutine  LOCNOD. 

For  purposes  of  efficiency,  the  cell  structure  should  extend  only  over 
the  domain  in  which  interaction  is  expected.  This  includes  the  domain  occu¬ 
pied  by  the  master  elements  in  the  undeformed  configuration  and  any  part  of 
the  domain  into  which  master  elements  are  anticipated  to  move  as  a  result  of 
the  interaction.  Although  the  optimal  relation  between  cell  size  and  master 
element  size  has  to  be  determined,  a  cell  should  span  at  least  two  element 
lengths  in  each  direction  to  reduce  the  number  of  elements  which  occupy  more 
than  one  cell.  A  cell  structure  consisting  of  3  x  3  x  1  cells  in  the  x,  y  and 
z  directions  as  shown  in  Fig.  3  has  been  found  quite  effective  in  the 
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identify  all  slave  nodes  in  a  cell,  and  then  in  treating  an  element  to 
identify  the  cell  number  in  which  it  is  located  and  to  check  the  slave  nodes 
which  occupy  the  same  cells  in  the  interaction  algorithm. 


The  stepe  of  the  interaction  procedures  within  the  structure  of  the 
complete  algorithm  is  shown  in  Table  4.  As  can  be  seen,  the  locations  of  the 
slave  nodes  are  determined  during  the  time  integration  of  the  velocities. 
Penetration  of  master  elements  by  slave  nodes  is  then  determined  in  a  loop 
over  all  master  elements.  This  loop  must  precede  the  element  loop  in  which 
new  nodal  forces  are  determined  because  the  interaction  algorithm  modifies  the 
velocities  of  the  slave  nodes  and  master  element  nodes. 


Table  4 

1 .  initial  conditions:  velocities  and  positions  of  all  nodes 

2.  integrate  velocities  to  obtain  new  displacements 

3. *  determine  the  cell  locations  of  all  slave  nodes 

4. *  for  each  master  element: 

4.  compute  surface  normal  vectors  and  assemble  into  global  array 

4b.  determine  cells  in  which  element  is  located 

4c.  by  checking  all  slave  nodes  in  these  cells,  determine  if  any  slave 
nodes  are  in  the  element 

4d.  if  a  slave  node  is  in  the  element,  move  it  back  to  an  outside  surface 
and  transfer  the  momentum  to  the  element  nodes,  which  modifies  its 
nodal  velocities 

5.  for  each  element: 

5a.  compute  strain- rates  from  the  nodal  velocities  and  stress- rates  from 
the  constitutive  equations 

5b.  integrate  stress- rates  to  obtain  new  stresses  and  compute  nodal 
forces 

5c.  assemble  nodal  forces  into  global  array 

6.  find  nodal  accelerations  from  equations  of  motion 

7.  integrate  accelerations  to  find  new  velocities;  go  to  2 

*  denotes  steps  which  pertain  to  the  interaction  algorithm 

2.2  Cell  Structure  and  Location  of  Nodes 


The  cell  structure  which  is  used  to  identify  the  slave  nodes  and  master 
elements  for  which  interaction  is  possible  is  shown  in  Fig.  3.  The  cells  are 
uniform  in  size  in  all  directions  and  their  edges  are  parallel  to  the  x,  y  and 
z  coordinates. 

The  cell  domain  is  described  by  the  following: 

NCX>  (2.2.1) 


Section  2 


INTERACTION  ALGORITHMS 


2.1  Overview 

The  basic  purpose  of  this  algorithm  is  to  treat  the  interaction  of  two 
bodies  with  eroding  elements.  Eroding  elements  are  elements  which  are 
destroyed  during  the  course  of  the  computation  because  of  very  high  strains, 
which  implies  that  they  represent  elements  of  the  target  or  penetrator  which 
have  ceased  to  play  a  significant  role  in  the  physics  of  the  problem. 

Algorithms  with  eroding  elements  in  three  dimensions  are  very  difficult 
to  treat  within  the  conventional  framework  of  master  and  slave  contact  sur¬ 
faces.  It  is  difficult  to  design  algorithms  to  redefine  the  contact  surfaces 
when  elements  are  destroyed,  particularly  in  three  dimensions.  Moreover,  the 
corners  which  are  created  in  the  surfaces  by  the  erosion  of  elements  present 
severe  algorithmic  difficulties. 

Therefore,  this  new  algorithm  uses  a  concept  of  slave  nodes  and  master 
elements.  One  of  the  two  bodies,  usually  the  projectile,  is  defined  by  the 
nodes,  hereafter  called  slave  nodes;  the  second  body  is  defined  by  elements, 
called  master  elements.  If  a  slave  node  is  found  to  be  inside  a  master 
element,  it  is  then  brought  back  to  an  outside  surface.  In  order  to  determine 
the  outside  surfaces,  a  set  of  normal  vectors  is  assembled  for  all  master 
elements  as  described  in  Section  2.4.  Because  of  the  character  of  the 
assembly  process,  a  nonzero  normal  vector  will  only  result  on  outside  surfaces 
and  will  provide  an  effective  average  normal  to  the  surface. 

The  mechanics  of  the  interaction  of  the  two  bodies  is  executed  completely 
through  the  interaction  of  the  slave  nodes  with  the  master  elements.  The 
rules  of  this  interaction  are  as  follows: 

i)  Slave  nodes  are  not  permitted  to  penetrate  master  elements. 

ii)  whenever  penetration  of  a  slave  node  into  a  master  element  is  detected, 
the  slave  node  is  returned  to  the  surface  of  the  element  it  has  penetrated 
and  the  associated  loss  of  momentum  is  transferred  to  the  appropriate 
nodes  of  the  master  element.  If  a  check  on  nodal  normals  shows  that  this 
is  not  an  exterior  surface,  the  node  is  moved  to  the  appropriate  edge. 

The  efficacy  of  this  procedure  hinges  strongly  on  the  use  of  explicit  time 
integration,  since  stability  requirements  then  limit  the  time  step  so  that  the 
master  element  which  is  penetrated  effectively  represents  the  zone  of  the 
interaction.  Because  of  the  large  number  of  slave  nodes  and  master  elements 
involved  in  this  process,  special  techniques  are  needed  to  quickly  identify 
the  slave  nodes  and  master  elements  between  which  interaction  is  possible. 

This  is  accomplished  by  using  a  cell  structure  which  is  fixed  in  space.  Cells 
are  considered  to  be  substantially  larger  than  elements  and  so  may  include 
many  master  elements  and  slave  nodes.  The  basis  of  the  procedure  is  to 
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At  -  ■§  [d  +  52)^  -  c]  (1.6.2) 

u 

where  5  is  the  fraction  of  critical  damping  in  the  maximum  frequency.  This 
time  step  is  always  smaller  than  the  stable  time  step  for  linear  problems. 
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(1.5.2) 

In  the  above, 

the  vectors  h 
-a 

are  defined  in 

Table 

3. 

Table  3 

.  The 

h  and  s 
** 

vectors  for 

the  hexahedron 

I  - 

1 

2 

3 

4 

5 

6 

7 

8 

h1I 

+1 

+  1 

-1 

-1 

-1 

-1 

+1 

+  1 

h2I 

+  1 

-1 

-1 

+1 

-1 

+1 

+1 

-1 

h3I 

+1 

-1 

+  1 

-1 

+  1 

-1 

+  1 

-1 

h4I 

-1 

+1 

-1 

+1 

+  1 

-1 

+  1 

-1 

A1I 

-1 

+1 

+  1 

-1 

-1 

+1 

♦  1 

-1 

A2I 

-1 

-1 

+1 

+1 

-1 

-1 

♦  1 

+1 

A3I 

-1 

-1 

-1 

-1 

+  1 

+1 

+1 

+1 

The  hourglass  strain- rates  give  the  rates  of  the  generalized  hourglass 
stresses 


*io 


°i)  qi« 


(1.5.3) 


where  the  rate  of  Qia  is  related  to  its  frame- invariant  rate  of  Q*a  by 


2ia  “  2 


ia  +  Wij 


(1.5.4) 


(1.5.5) 


The  constants  are  given  by 
Gij  ■  *  1  sij 

where  w  is  an  estimate  of  the  maximum  frequency  and  e  a  user- control led 
parameter.  Values  of  0.03  to  0.10  are  recommended  for  e. 


1 .6  Element  Frequency  and  Stable  Time  Step 

The  stable  time  step  for  the  element  is  computed  by  using  the  following 
upper  bound  for  the  maximum  element  frequency  taken  from  Ref.  [4] 


u2  <  u)2 
max 


^  Bi!  Bil 


(1 .6.1 ) 


where  c  is  the  dilatational  wavespeed.  The  time  step  computed  by 


The  advantage  of  this  form  over  Bq.  (3.1)  is  apparent  in  Bq.  (3.6),  which 
shows  that  the  last  four  terms  of  B.  need  not  be  computed.  On  the  other 
hand,  Bq.  (3.1)  applies  to  degenerate  hexahedra  such  as  tetrahedra,  whereas 
Eq.  (3.5)  does  not. 


1.4  Stress- Strain  Relationship 


The  element  can  use  any  nonlinear  stress- strain  law,  but  for  an  aniso¬ 
tropic  stress- stirain  law  the  material  matrices  must  be  updated  as  indicated  in 
Ref.  [5],  which  is  indicated  by  Bq.  (1.4.4). The  measure  of  stress  in  this 
element  is  the  Cauchy  stress,  or  physical  stress,  and  the  measure  of  deform¬ 
ation  is  the  velocity  strain,  also  known  as  'rate-of-deformation.  This  will 
usually  be  called  the  strain- rate.  In  order  to  maintain  frame- invariance  in  a 
stress- strain  law  based  on  these  measures,  a  frame- invariant  stress  rate  must 
be  used.  Here,  the  Jaumann  rate  is  used.  If  the  material  stress  rate  is 
related  to  the  strain  rate  by 


Cijki 


(1.4.1) 


then  the  rate  of  change  of  the  stress  is  given  by 


+  Wik  °kj 


+  Wjk  °ki 


~ij  ij 

For  any  isotropic  material,  Bq.  (4.1)  can  be  written  as 


(1 .4.2) 


^k 


+  2v  e 


ij 


(1.4.3) 


where  X  and  |i  are  the  Lame  constants. 


This  form  of  the  stress- strain  relations  is  frame- invariant  only  if  the 
material  is  and  remains  isotropic.  For  anisotropic  materials,  the  C  matrix 
must  be  updated  as  follows  [5] 

*  7 

c. .. ,  ■  c....  +  w.  c  ...  +  w..  c.. , .  +  w.  c..  .  +  w  .  c...  . 

ijki  ij)ci  ia  ajki  ]b  ibki  kc  ijcl  Jtd  ijkd  . 


1 .5  Hourglass  Control 


(1  .4.4) 


When  the  nodal  forces  of  the  hexahedron  are  evaluated  by  one-point 
quadrature,  it  possesses  12  spurious  zero-energy  modes  or  hourglass  modes. 

The  modes  occur  independently  in  the  x,  y  and  z  directions.  The  four  modes  in 
the  x-direction  are  shown  in  Fig.  2;  the  modes  in  the  y  and  z  directions  can 
be  envisioned  by  simply  replacing  the  x-axis  by  the  y  or  z  axes. 


To  control  these  modes,  12  additional  generalized  strain- rates  are 
defined  by 

«ia  "  Yal  vil  (K5-1> 

where  the  range  of  Greek  letters  is  4,  and 
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In  the  centroidal  form,  the  B  matrix  is  evaluated  at  the  centroid  of  the 
element,  that  is,  the  point  £  -  n  »  C  -  0.  In  this  form,  the  B  matrix  is 
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(1.3.5) 

where  J  is  the  Jacobian  at  the  centroid.  The  terms  B  _  to  B  .  are  then 
obtained  by  x5  x8 


(1.3.6) 


The  terms  of  g  associated  with  y  and  z,  i.e.  the  second  and  third  rows,  are 
obtained  by  the  same  permutations  of  x,  y,  and  z  on  the  right  hand  side  as 
indicated  in  Table  2. 


The  Jacobian  is  given  by 
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and  A.,  are  given  in  Table  3. 
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(1 .3.8a) 
(1.3.8b) 


(1  .3.8c) 


Table  1 , 

,  Permutations  of  node  numbers  for  generating 
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+ 

y7(Z68  + 

Z24)  + 

y8  Z74  +  y6  Z27 

• 

The 

other  three 

rows  are  obtained  by  interchanging  x,  ; 

y  f  z 

according  to 

Table  2. 

Table  2.  Permutations  of  coordinates  for 

generating 

rows  2  and  3  from  row 

1 

Row 

1 

y 

z 

2 

z 

X 

3 

X 

y 

For 

example 

12V  B33 

-  ( 

X4  (y81  + 

y72)  + 

X1  y42  +  X2(yi6 

+  y47 

) 

4 

'  X7(y68  + 

y24)  + 

X8  y74  +  X6  y27 

]. 

(1 .3.3) 

The 

volume  of  the  element 

is  obtained  by 
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contour  are  considered. 


J?=x42xx31  (2.4.1) 

where  (x)  designates  a  vector  cross-product.  This  vector  n  is  normalized  and 
then  assembled  into  the  global  arrays  of  the  normals  of  4  nodes  by  adding  each 
component  of  the  vector  p  to  the  existing  vector  in  the  nodal  array.  When  the 
contribution  of  each  element  has  been  added  to  the  nodal  arrays,  the  procedure 
is  complete.  The  procedure  can  be  summarized  as  follows: 

1 .  Initialization  at  beginning  of  time  step:  zero  the  vectors  p  for  all 
nodes  I . 

2.  For  each  element  e  in  the  set  of  master  elements: 

2.1  For  each  side  of  the  element,  compute  p  ? 

2.2  add  n  into  the  arrays  n  of  the  four  nodes  which  define  the  side. 


Remarks: 

1 .  In  order  to  avoid  difficulties  with  master  elements  at  the  edge  of  the 
master  element  domain,  the  master  elements  should  be  defined,  to  extend 
sufficiently  far  so  that  the  edge  master  elements  are  not  penetrated  from 
the  sides. 

2.  For  nodes  of  master  elements  which  lie  on  a  plane  of  symmetry,  the 
components  of  the  normal  vector  which  do  not  lie  in  the  plane  of  symmetry 
are  set  to  zero. 

3.  For  master  elements  which  lie  on  the  bottom  of  the  target  (opposite  to  the 
original  position  of  the  projectile),  the  normals  to  the  bottom  surface 
are  omitted  to  avoid  driving  the  slave  nodes  in  the  wrong  direction, 

2 .5  Adjustment  of  Slave  Node  Positions 

Once  it  is  determined  that  a  slave  node  has  penetrated  a  master  element, 
it  is  necessary  to  adjust  its  position  consistent  with  the  fact  that  its 
normal  momentum  has  been  transferred  to  the  target.  This  procedure  involves 
two  steps: 

1 .  using  the  normal  vectors  associated  with  the  element,  determine  in  which 
direction  the  position  of  the  slave  node  has  to  be  adjusted? 

2.  displace  the  slave  node  to  an  outside  surface  in  the  direction  of  the 
normal  p? 

3.  if  the  surface  to  which  a  slave  node  is  brought  is  not  an  outside  surface, 
bring  the  slave  node  back  to  an  edge  of  the  surface. 

The  procedure  is  implemented  as  follows.  The  average  normal  of  the 
element  is  found  by 

8 

8  -  (  I  ffJ  /  1  1  (2.5.1 ) 

1-1 

where  the  division  by  "H  I"  designates  normalization  of  the  vector,  the  norm 
is  defined  in  Eq.  (3. 3). Let  the  current  coordinates  of  the  slave  node  be  x  . 
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Then  the  node  is  displaced  by  the  procedure 


~n  3  ~o  +  n  ~  (2‘5‘2) 

where  n  is  an  undetermined  parameter  n  >  0. 

The  magnitude  of  n  is  determined  by  checking  which  of  the  sides  of  the 
hexahedron  is  intersected  by  the  line  of  Bq.  (2.5.2).  This  is  accomplished  as 
follows.  Each  side  is  subdivided  into  2  triangular  surfaces  (note  that  this 
is  only  an  approximation  to  the  surface  of  an  isoparametric  hexahedron  but  it 
simplifies  computations  enormously).  By  taking  3  nodes  of  the  surface  in 
turn,  the  surface  is  defined  by 

Xi“XiICX  1-1  to  3  (2.5.3) 

+  52  +  C3  -  1  .  (2.5.4) 

The  intersection  of  the  line  defined  by  Eq.  (5.2)  and  the  plane  is 
determined  by  solving  the  4  equations  in  4  unknowns  represented  by  Eqs.  (5.3) 
to  (5.4).  The  solution  is  given  by 
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A  particular  triangular  surface  is  intersected  by  the  parametric  ray  of 
Eq.  (5.2)  if  and  only  if 


h  >  0 


(2.5.7) ' 


0  <  Cj  <  1  for  I  *  1  to  3  •  (2.5.8) 

Once  the  surface  on  which  the  slave  node  is  projected  is  determined,  the 
surface  is  checked  to  ascertain  whether  it  is  an  outside  surface.  This  is 
done  by  checking  whether  the  4  normals  of  the  nodes  of  the  surface  are  non 
zero.  If  this  check  fails,  the  node  is  projected  to  an  edge  of  the  surface  as 
shown  in  Fig.  7. 


The  equations  which  govern  this  realignment  of  the  3lave  node  are  the 

following.  Let  x  be  given  by 
~n 

x  »  x  +  At  vold  (2.5.9) 

~n  ~o  ~ 

where  xq  is  the  position  of  the  slave  node  at  the  beginning  of  the  time  step 
as  shown  in  Fig.  7.  The  edges  of  the  side  are  then  considered  in  turn  and  its 
nodes  are  generically  identified  as  1  and  2,  with  the  vector  connecting  them 
denoted  by  .  The  previously  computed  new  position  of  the  node  is  denoted 
by  x  where 

x  »x  +  A  r  .  (2.5.10) 

~m  ~n  ~ 


The  node  is  then  repositioned  on  the  intersection  of  the  line  with  plane 

defined  by  the  vectors  x  and  x  .  The  equations  to  be  solved  are 

~m  ~n 


~10 


x,,  £  -  x  -  x 

~21  1  ~on  2  ~i 


mn 


(2.5.11  ) 


where  5^,  i  *  1  to  3,  are  the  unknowns.  If  the  edge  $2i  is  the  correct 
one  £  must  satisfy  0  <  £1  <1.  Otherwise,  another  edge  of  the  surface  is 
considered. 


The  final  position  is  denoted  by  x^  which  is  given  by 

"*■  ^21  (2. 5.1 2a) 

where  £1  is  the  above  solution.  The  reposition  vector  is  then  redefined  by 


Ar  »  x. 


x,  -  x 


(2.5.12b) 


Procedure  for  repositioning  a  slave  node  when  the  normal  projection 
leaves  it  inside  the  target;  node  B  illustrates  a  case  where  velocity 
is  increased  by  original  repositioning. 


If  a  slave  node  is  on  a  plane  of  symmei 
to  the  plane  of  symmetry  is  now  set  to  zero 


the  component  normal  of  Ar 


One  situation  which  can  foil  even  this  check  and  allow  a  slave  node  to 
remain  within  the  target  is  shown  in  Fig.  8.  Here  the  slave  node  has 
penetrated  element  with  nodes  numbered  1  to  8.  The  repositioning  will  place 
the  node  on  the  surface  defined  by  nodes  1  to  4,  and  because  a  single  intact 
element,  defined  by  nodes  1  to  4  and  9  to  12  lies  above  the  intact  plane,  the 
normals  will  be  nonzero  on  each  of  the  four  nodes,  1  to  4.  Therefore,  the 
node  will  remain  mi  this  surface.  However,  if  all  projectile  nodes  are  slave 
nodes,  a  subsequent  slave  node  should  engage  the  element  defined  by  nodes  1  to 
4  and  9  to  12. 


Once  the  new  position  of  the  slave  node  is  determined  the  change  in  its 
velocity  is  computed  by 


Av  ■  Ar/At 


(2.5.13a) 


where  Ar  is  the  total  displacement  of  the  slave  node  needed, 
the  slave  node  is  then  modified  by 


new 


old 


+  Av 


The  velocity  of 


(2.5.13b) 


In  unusual  circumstances  it  is  possible  for  the  normal  of  an  element  to 
form  an  angle  of  more  than  90°  with  the  slave  node  velocity.  This  situation 
is  illustrated  in  Fig.  8,  for  node  B,  where  a  is  the  angle  between  the 
velocity  vector  and  the  normal.  If  the  previously  described  procedures  are 
used,  the  repositioning  of  the  slave  node  will  increase  its  velocity  and  its 
kinetic  energy  as  shown;  note  the  new  velocity  vector  is  longer  than  the 
original  velocity  vector.  This  repositioning  is  of  course  completely  contrary 
to  physical  laws,  since  it  increases  the  kinetic  energy  of  the  system. 

new 

In  order  to  avoi^this,  the  new  velocity  v  is  compared  with  the 
original  velocity,  v  .  It  is  required  that  the  following  inequality  be 
satisfied  by  these  velocities 

I  v  I  <  I  v  I  *  (2.5.14) 

new 

If  this  is  not  satisfied,  the  node  is  moved  back  along  the  vector  v  until 
its  velocity  satisfies  the  following  equation 


1  vn*W  I2  +  I  A  r/At  I2  -  »  vold  I2  •  (2.5.15) 

new 

This  implies  that  the  angle  between  A  r  and  v  will  be  a  right  angle. 
Although  this  procedure  will  usually  position  the  slave  node  outside  of  the 
target,  it  insures  that  energy  is  not  generated  by  the  procedure.  In 
subsequent  time  steps  the  slave  node  will  again  penetrate  a  master  element  so 
the  procedure  is  not  harmful. 


The  momentum  loss  associated  with  this  adjustment  is  MAv,  where  M  is  the 
mass  of  the  slave  node.  This  mass  is  now  transferred  to  the  nodes  in  contact 
with  the  2  triangles  on  the  penetrated  side.  The  formula  used  is 


*  H  T  - 

Av  »  — —  n  n  Av  (no  sum  on  J) 

J  TM  ^  J  ^  ~ 

J 

where 

r  T  - 
r  -  l  Bt  P 
1-1 


(2.5.16) 


(2.5.17) 


This  formula  apportions  the  momentum  to  the  nodes  according  to  how  strongly 
their  vectors  point  in  the  direction  of  the  interface  normal  n. 
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Section  3 


NUMERICAL  EXAMPLES 


The  sample  problems  include  a  very  simple,  small  problem  which  can  be 
used  to  quickly  check  the  performance  of  the  program  in  a  new  installation  and 
a  moderate  size  problem  and  large  scale  problem.  In  each  case,  a  table  of  key 
problem  parameters,  a  printout  of  the  card  images  of  the  input  data  deck  and 
some  computer  graphics  of  the  evolution  of  the  response  are  given.  The 
execution  speed  of  the  program  is  100  elements  per  time  step  per  CPU  second  on 
a  CYBER  170/730,  1000  elements  per  time  step  on  a  CDC  7600. 

Example  1 

The  first  example  involves  a  simple  copper  projectile  consisting  of  24 
elements  striking  a  target  modeled  by  48  elements.  The  material  properties 
and  dimensions  are  given  in  Table  5.  The  evolution  of  the  problem  is  shown  in 
Fig.  9.  As  can  be  seen  from  Fig.  9,  because  of  the  small  size  of  the  target, 
although  erosion  commences,  subsequent  momentum  transfer  causes  the  target  to 
move  away  from  the  projectile. 

Example  2 


The  second  example  is  a  problem  of  moderate  scale,  involving  88  elements 
in  the  projectile  and  500  elements  in  the  target.  Table  6  gives  the  problem 
parameters  and  card  images  of  the  data.  The  evolution  of  the  problem  is  shown 
in  Fig.  10.  Erosion  is  only  evident  in  the  projectile  for  the  first  77  time 
steps  (cycles)  which  were  run. 

Example  3 

Example  3  is  taken  from  Ref.  [7].  It  consists  of  a  copper  rod  striking  a 
steel  plate  at  2000  meter/second.  The  projectile  is  modeled  with  414 
elements,  the  target  with  1014  elements.  The  copper  is  assigned  a  failure 
strain  of  i  =  2.0.  A  complete  listing  of  material  and  problem  parameters  is 
given  in  TaEle  7.  Note  that  arbitrary  erosion  can  occur  in  both  the  target 
and  projectile. 

The  simulation  is  shown  in  Fig.  11.  Note  that  the  projectile  starts 
jetting  in  the  positive  x-direction  early  in  the  simulation.  These  large 
3hears  result  in  rapid  erosion  of  the  projectile.  Subsequently,  large 
deformations  in  the  target  result  in  shear  failure  in  the  target.  Erosion 
takes  place  in  both  the  target  and  projectile. 

In  Fig.  11,  gaps  often  appear  to  develop  between  the  target  and 
projectile.  This  is  partially  a  result  of  the  use  of  a  two  dimensional  plot 
of  the  plane  of  symmetry,  which  cannot  show  the  contact  between  the  projectile 
and  target  away  from  the  plane  of  symmetry. 


TABLE  5 


Projectile 

shape 

dimensions 
density 
bulk  modulus 
shear  modulus 
yield  stress 
ultimate  stress 
initial  velocity 

Target 

shape 

dimensions 
density 
bulk  modulus 
shear  modulus 
yield  stress 
ultimate  stress 


Parameters  and  input  for  Example  1 


:  rod 

:  3  in  long,  radius  0.6  in 

:  0.000831  lb-sec2/in4 

:  20,739,000  psi 

s  6,380,000  psi 

s  20,300  psi 

:  65,300  psi 

s  x- component  -  2588.0  in/sec 

z-component  -  9659.0  in/sec 


:  plate 

:  6  in  x  3  in  x  0.5  in  (half  plate) 

s  0.0005  lb-sec2/in4 
:  24,200,000  psi 

:  9,300,00  psi 

:  160,000  psi 

:  160,100  psi 


initial  velocity 


0 


Table  5  Continued 


l  1  1 

I  1  PROJECTILE  MATERIAL 
.00073  410000.0  ?3*lO00O. 


30000. 

250000. 

310000. 

0.3 

qqqqqqq. 

0.0 

1.0 

999. 

090. 

2.0 

0.0 

>  l 

TARGET  MATERIAL 

.00073 

410000.0 

2773  0000 . 

30000. 

160000. 

185000. 

0.3 

qqqqqqq. 

1.0 

999. 

099  . 

0.1 

0.0 

1.0 

1.0 

1  .0 

0.0 

0.2 

15.0 

0.0 

9  0 

4 

0 

3.0 

0.0 

1.0 

0.6 

O.o 

0.6 

0.0 

1.0 

1.0 

1.0 

4  1 

?  4 

3 

3 

4 

1 

1.2 

0.3 

1.2 

1.0 

1 

-3.0 

0.0 

0.0 

3.0 

3.0 

-0.5 

?  1 

3 

1 

0 

1 

0 

V  1 

3  3 

2 

61 

2 

1 

1 

8  3 

1 

l  3 

1 

61 

1 

60 

60 

0 

0.0 

1.0 

1 

5  4 

1 

9 

1 

373  2,5 

0.0 

-62633.1 

0.000001 
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Table  8 


MAIM  ROUTINE  INPUT  DATA 


DESCRIPTION  CARD  (12A6) 
j Descrip uion  of  Problem 


IDENTIFICATION  CARD  (315 ,5X,F10 .0,E1 5 .8, 5X,F10 .0 ) 


CASE  |  CYCLE 

IPRES 

1 - 1 

CPMAX 

EMAX  I 

HGMAX  { 

INTEGRATION  TINE  INCREMENT  CARD  AND  CELL  STRUCTURE  (7F10.0) 
DTHAX  1  DTMIN  1  SSF  I  TMAX  I  XDIS  I  YD IS  j  ZD IS  I 


CELL  PARAMETER  CARD  (6E10.3) 


| SXMAX  l  SXMIN 

SYMAX 

SYMIN 

SZMAX 

SZMIN 


If  NX1  *  1,  NX2  -  HLX,  NY 2  =*  NLY  and  NZ2  -  NLZ-1 ,  all  elements  in  the  target 
are  master  elements. 

Main 


Identification 

Card  (315,  5X,  F10.0,  E15 

.8,  FI 0.3)  -  Add  one  additional 

,v 

parameter  SFAIL  in  column  51-60. 

SFAIL:  hourglass  failure  criterion; 

(10.0~12.0  is  recommended) 

Integration  Time  Increment  Card  (8F10.0) 

- 

A-l 

< 

column 

variable  name 

description 

1  -  10 

DTMAX 

maximum  time  step 

11-20 
21  -  30 

DTMIN 

SSF 

minimum  time  step 
time  step  safety  factor, 

<  1 .0 

-  « 

31  -  40 

TMAX 

time  problem  is  allowed  to 

i  run 

41  -  50 

HRC0N 

hourglass  control  factor, 

0.05  to 

51  -  60 

XDIS 

0.2  is  recommended, 
size  in  x-direction  for  a 

cell 

61  -  70 

YD  IS 

size  in  y-direction  for  a 

cell 

•„ 

71  -  80 

2D  IS 

size  in  z- direction  for  a 

cell 

hr* 

Cell  Parameter 

Card  (6E10.3)  -  See  Section  2.2 

*  _  j 

This  card 

should  follow  right  after 

the  Integration  Time  Increment  Card. 

wmJm 

I*  ^ 

column 

variable  name 

description 

1  -  10 

SXMAX 

x_ _ 

\~4 

11-20 

SXMIN 

max 

X  i 

21  -  30 

31  -  40 

SYMAX 

SYMIN 

min 

ymax 

ymin 

41  -  50 

SZMAX 

z 

51  -  60 

SZMIN 

max 

z  , 

•/•I 

min 

Only  the  metallic  material  (material  code  -  1)  can  be  used.  The  input  format 


for  the  Main  Program  is  summarized  in  Table  8. 


Postprocessor 

Only  the  plane  of  symmetry  can  be  plotted  when  the  hexahedra  and  erosion 
features  are  used. 
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Section  4 


INPUT  FORMAT 


The  input  format  is  almost  identical  to  that  of  the  original  EPIC- 3  code 
[6,7],  The  major  differences  are  that  the  cell  description  has  been  added  to 
the  integration  time  increment  card  and  a  subsequent  card,  and  there  are  some 
restrictions  on  the  features  which  cam  be  used  with  the  hexahedron  and 
erosion . 

The  following  restrictions  apply  to  the  hexahedron  and  erosion  interface: 

1.  the  erosion  interface  can  not  be  used  with  tetrahedral  elements; 

2.  the  anisotropic  material  cannot  be  used  with  the  hexahedra. 

Preprocessor 

The  master  elements  must  be  identified  through  the  element  description 
cards.  This  is  accomplished  by  specifying  MIDEN  in  columns  61-65  for  the 
composite-element-description  cards  or  columns  31-35  for  any  other  element 
description  cards.  MIDEN  is  specified  as  follows: 

0  (or  blank)  if  the  element  is  not  a  master  element 

MIDEN  =*  1  if  the  element  is  a  master  but  not  on  the  bottom  layer 

2  if  the  element  is  a  master  on  the  bottom  layer  of  the  target 

For  a  plate  target,  MIDEN  can  be  generated  automatically  if  the  master 
elements  are  to  occupy  a  regular  domain  in  the  target.  Recall  that  the  flat- 
plate  target  consists  of  NLX  layers  of  elements  in  the  x-direction,  NLY  layers 
in  the  y-direction,  and  NLZ  layers  in  the  z-direction.  We  assign  numbers  to 
the  layers  beginning  with  lowest  value  of  the  coordinate  and  in  the  direction 
in  which  the  coordinate  increases.  The  master  elements  identification  (MIDEN) 
can  then  be  generated  automatically  whenever  all  elements  between  layers  NX1 
and  NX2  in  the  x-direction,  between  layers  1  and  NY2  in  the  y-direction,  and 
layers  1  and  NZ2  in  the  z-direction  are  master  elements.  In  that  case, 
columns  31-35  of  the  flat-plate  element  description  card  are  as  follows: 

Flat  Plate  Description  Card 


Columns 

Variable  Name 

Description 

31-35 

MIDEN 

leave  blank  when  automatic  generation 
of  MIDEN  is  to  be  used. 

36-40 

NX1 

lowest  layer  number  for  master  elements 
in  x-direction 

41-45 

NX  2 

largest  layer  number  for  master 
elements  in  x-direction 

46-50 

NY2 

largest  layer  number  for  master 
elements  in  y-direction 

51-55 

NZ2 

largest  layer  number  for  master 
elements  in  z-direction 

TIME  =.000015092 
CYCLE  =62 
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Fig.  11a  Evolution  of  mesh  on  plane  of 
symmetry  for  example  3. 
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Table  7  Continued 
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TABLE  7 


Parameters  and  input  for  Example  3 


Projectile 

shape 

• 

rod  with  a  round  nose 

dimensions 

: 

4.9  in  long,  0.5  in  radius 

density 

: 

0.000831  lb-sec2/in4 

bulk  modulus 

: 

20,739,000  psi 

shear  modulus 

• 

6,380,000  psi 

yield  stress 

: 

20,300  psi 

ultimate  stress 

s 

65,300  psi 

initial  velocity 

•* 

x-component  55660.0  in/sec 
z-component  -55660.0  in/sec 

Target 

shape 

dimensions 
density 
bulk  modulus 
shear  modulus 
yield  stress 
ultimate  stress 
initial  velocity 


plate 

7.9  in  x  3.95  in  x  0.375  in  (half  plate) 
0.000734  lb-sec2/in4 
24,200,000  psi 
9,300,000  psi 
160,000  psi 
203,000  psi 
0 


45 


TIME  =.000012519 
CYCLE  =77 


TIME  =.000007129 


Fig.  10b  Evolution  of  mesh  for  example 


Table  6  Continued 
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TABLE  6 


Projectile 

shape 

dimensions 
density 
bulk  modulus 
shear  modulus 
yield  stress 
ultimate  stress 
initial  velocity 

Target 

shape 

dimensions 
density 
bulk  modulus 
shear  modulus 
yield  stress 
ultimate  stress 


Parameters  and  input  for  Example  2 


:  rod 

;  4.04  in  long,  0.201  in  radius 

:  0.00073  lb-sec2/in4 

:  23,810,000  psi 

s  1 1 ,630,000  psi 
:  250,000  psi 

:  310,000  psi 

:  x-component  56155.0  in/sec 
z-component  -32421.0  in/sec 


s  plate 

:  5.6  in  x  0.7  in  x  1  in  (half  plate) 

s  0.00073  lb-sec2/in4 
:  27,780,000  psi 

:  1 1 ,360,000  psi 

:  160,000  psi 

:  185,000  psi 


initial  velocity 


0 
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This  Laboratory  undertakes  a  continuing  effort  to  improve  the  quality  of  the 
reports  it  publishes.  Your  conunents/answers  to  the  items/questions  below  will 
aid  us  in  our  efforts. 

1.  BRL  Report  Number _ Date  of  Report _ 

2.  Date  Report  Received _ 

3.  Does  this  report  satisfy  a  need?  (Comment  on  purpose,  related  project,  or 

other  area  of  interest  for  which  the  report  will  be  used.) _ 


4.  How  specifically,  is  the  report  being  used?  (Information  source,  design 
data,  procedure,  source  of  ideas,  etc.) _ 


5.  Has  the  information  in  this  report  led  to  any  quantitative  savings  as  far 
as  man-hours  or  dollars  saved,  operating  costs  avoided  or  efficiencies  achieved, 
etc?  If  so,  please  elaborate. _ 


6.  General  Comments.  What  do  you  think  should  be  changed  to  improve  future 
reports?  (Indicate  changes  to  organization,  technical  content,  format,  etc.) 
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Address 


City,  State,  Zip 

7.  If  indicating  a  Change  of  Address  or  Address  Correction,  please  provide  the 
New  or  Correct  Address  in  Block  6  above  and  the  Old  or  Incorrect  address  below. 
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